A trivial SOS decomposition example

Contributed by: votroto Adapted from: (Blekherman et al., 2012; Examples 3.25)

using DynamicPolynomials
using SumOfSquares
import CSDP

The polynomial p = x^2 - x*y^2 + y^4 + 1 is SOS. We can, for example, decompose it as p = 3/4*(x - y^2)^2 + 1/4*(x + y)^2 + 1, which clearly proves that p is SOS, and there are infinitely many other ways to decompose p into sums of squares.

We can use SumOfSquares.jl to find such decompositions.

First, setup the polynomial of interest.

@polyvar x y
p = x^2 - x*y^2 + y^4 + 1

\[ 1 + x^{2} - xy^{2} + y^{4} \]

Secondly, constrain the polynomial to be nonnegative. SumOfSquares.jl transparently reinterprets polyonmial nonnegativity as the appropriate SOS certificate for polynomials nonnegative on semialgebraic sets.

model = SOSModel(CSDP.Optimizer)
@constraint(model, cref, p >= 0)

\[ (1) \cdot 1 + (1) \cdot x^{2} + (-1) \cdot xy^{2} + (1) \cdot y^{4} \text{ is SOS} \]

Thirdly, optimize the feasibility problem!

optimize!(model)

CSDP 6.2.0
This is a pure primal feasibility problem.
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 9.00e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  2.7659722e+01
Iter:  2 Ap: 1.00e+00 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  1.2562719e+01

Lastly, recover a SOS decomposition. In general, SOS decompositions are not unique!

sos_dec = sos_decomposition(cref, 1e-4)

(-0.7423566403018642·1 - 0.5916342145954683·x + 0.9492757372229985·y^2)^2 + (1.7427686231333871e-16·1 + 1.1201589623631873·y + 1.5284610719333535e-16·x + 2.777860369511558e-16·y^2)^2 + (0.6232480104529257·1 - 0.782024243528594·x - 5.528165712762894e-17·y^2)^2 + (-0.24590350966628932·1 + 9.87415030380184e-17·y - 0.19597713808894557·x - 0.3144448675359998·y^2)^2

Converting, rounding, and simplifying - Huzza, Back where we began!

polynomial(sos_dec, Float32)

\[ 1.0 + 3.4187383 \cdot 10^{-16}y + 5.689893 \cdot 10^{-16}x - 8.881784 \cdot 10^{-16}y^{2} + 3.0372172 \cdot 10^{-16}xy + x^{2} + 5.602315 \cdot 10^{-16}y^{3} - xy^{2} + y^{4} \]

A deeper explanation and the unexplained 1e-4 parameter

p = x^2 - x*y^2 + y^4 + 1 can be represented in terms of its Gram matrix as

gram = gram_matrix(cref)

GramMatrix with row/column basis:
 SubBasis{Monomial}([1, y, x, y^2])
And entries in a 4×4 SymMatrix{Float64}:
  1.0                     0.0                …  -0.6273780504812865
  0.0                     1.254756100962573      0.0
  3.3723024372989114e-17  0.0                   -0.49999999999999994
 -0.6273780504812865      0.0                    1.0

keys_as_monomials(gram.basis)' * gram.Q * keys_as_monomials(gram.basis)

\[ 1.0 + 6.744604874597823 \cdot 10^{-17}x + x^{2} - 0.9999999999999999xy^{2} + y^{4} \]

where the matrix gram.Q is positive semidefinite, because p is SOS. If we could only get the decomposition gram.Q = V' * V, the SOS decomposition would simply be ||V * monomials||^2.

Unfortunately, we can not use Cholesky decomposition, since gram Q is only semidefinite, not definite. Hence, SumOfSquares.jl uses SVD decomposition instead and discards small singular values (in our case 1e-4).

This page was generated using Literate.jl.